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Electron transport in graphene is along the sheet but junction devices are often made by stacking different 
sheets together in a “side-contact” geometry which causes the current to flow perpendicular to the sheets within 
the device. Such geometry presents a challenge to first-principles transport methods. We solve this problem 
by implementing a plane-wave based multiple scattering theory for electron transport. This implementation im¬ 
proves the computational efficiency over the existing plane-wave transport code, scales better for parallelization 
over large number of nodes, and does not require the current direction to be along a lattice axis. As a first 
application, we calculate the tunneling current through a side-contact graphene junction formed by two separate 
graphene sheets with the edges overlapping each other. We find that transport properties of this junction de¬ 
pend strongly on the AA or AB stacking within the overlapping region as well as the vacuum gap between two 
graphene sheets. Such transport behaviors are explained in terms of carbon orbital orientation, hybridization, 
and delocalization as the geometry is varied. 


I. INTRODUCTION 

Graphene has some of the most fascinating electrical, ther¬ 
mal and mechanical properties 1 2 that promise to make it an 
important material for broad applications. However, to realize 
such a promise we need to learn more than its bulk properties. 
For example, most large-area graphene films are produced as 
polycrystalline sheets 3 6 containing multiple small domains, 
usually connected by one of two types of boundaries, end 
contacts where two domains connect within the same sheet 
with direct atomic bonding and usually referred to as grain 
boundaries 7 ®, and side contacts formed by stacking edge re¬ 
gions of the two graphene domains side by side with van der 
Waals force holding them together, which are also observed 
in recent experiment 12 . There is a great potential for useful 
devices 13 ® using side-contact junctions, in which the over¬ 
lapping region is the device region while the rest of the two 
graphene domains act as electrodes. First-principles trans¬ 
port study of either types of boundaries in graphene sheets 
can be challenging because of poor screening due to low di¬ 
mensionality. Moreover, the side-contact junctions present a 
particularly difficult problem because of its unaccommodating 
geometry for computational methods designed to deal with 
layer-structured systems. 

There are two basic approaches to apply the mesoscopic 
theory of Laudauer and Buttikei®® to study quantum trans¬ 
port of electrons within the first-principles method. The first 
approach is to use localized basis sets 17 20 which allow the 
calculation of the Green’s function of an electrode-device- 
electrode assembly and a straightforward transport calcula¬ 
tion based on the Green’s function method. However, local¬ 
ized basis sets do not work well for tunneling through large 
vacuum gaps that require a faithful description of vacuum 
electron wave functions. The second approach is based on 
the scattering theory of plane waves. It has been shown to 
be completely equivalent to the nonequilibrium Green’s func¬ 
tion method in the case of noninteracting electron ^ 21 l A rig¬ 
orous first-principles method 22 23 (Choi-Ihm method) based 


on scattering theory and pseudopotentials is implemented 
within the PWCOND part of the QUANTUM ESPRESSO 
package 2 ^ This code has become one of the standard tools 
for quantum transport studies 25 29 . Within this second ap¬ 
proach there is also a somewhat different implementation 
based on the Korringa-Kohn-Rostoker (KKR) band theory, 
alternatively called the multiple scattering theory, that is 
adapted to layer-structured systems and thus named the layer- 
KKR method 2130 . This method has been particularly suc¬ 
cessful in the study of spintronics ^ 1 33 . While localized basis 
is inadequate for the study of tunneling current between two 
graphene sheets in the side-contact junction, neither the Choi- 
Ihm method nor the layer-KKR method can be effectively ap¬ 
plied to this problem as well, as we will discuss below. This 
dissatisfaction compels us to search for a third implementation 
of the plane wave scattering method. 

Both the layer-KKR and Choi-Ihm methods use a two- 
dimensional plane-wave basis and divide the system under 
study into a stack of sufficiently thin slices along the transport 
direction. A generalized complex band structurd® or trans¬ 
mission matrices are then computed by stitching these slices 
together with appropriate boundary conditions. Each method 
has its own advantages and drawbacks. On the one hand, the 
layer-KKR method first solves the Green’s function of indi¬ 
vidual atomic layers and then stitches these layers together 
using a layer doubling technique based on multiple scattering 
theory. This approach is more efficient and yields the scatter¬ 
ing matrix for any part of the system of interest thus providing 
more information about the transport properties of the system. 
However, it has a serious drawback. It is implemented within 
the muffin-tin or the atomic sphere approximations (ASA), re¬ 
quiring that the space be divided into spheres around each 
atom within which the Kohn-Sham potential is spherically 
symmetric. On the other hand, the Choi-Ihm method does 
not require spherical approximations. It uses much thinner 
slices and stitches the slices together by matching boundary 
conditions of the wave functions between the slices. While 
this method is in principle rigorous, it is computationally ex¬ 
pensive. Moreover, it only yields transport information for 
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the whole system at the end without providing any scattering 
matrix of the individual parts of the system. 

Here we present a plane-wave multiple scattering algorithm 
that combines the advantages of the layer-KKR and Choi-Ihm 
methods, and show that it can be an efficient first-principles 
quantum transport method for tunnel junctions represented by 
the side-contact junction of graphene sheet discussed above 
and other devices characterized by low symmetry. We first 
completely reformulate the multiple scattering theory within a 
plane-wave basis, then implement an algorithm similar to the 
layer-KKR method but without the muffin-tin or ASA approx¬ 
imations. Additional computational speedup is realized by 
calculating the complex bands of each electrode using only its 
two-dimensional primitive unit cell and then folding the com¬ 
plex bands into the smaller Brillouin zone corresponding to 
the larger two-dimensional supercell of the scattering region 
(described in Section[Tl]), and by incorporating the method de¬ 
veloped by Srivastava et aim for low-symmetry nonorthog- 
onal lattices based on the three-dimensional Bloch theorem. 
While the impementation of both improvements is straightfor¬ 
ward, neither is available in the PWCOND code. As a first ap¬ 
plication, we use the new method to study conduction through 
a graphene/graphene side-contact junction. Our calculations 
show that there are significant differences in transmission co¬ 
efficients and their scaling with the overlapping area between 
AA and AB stacked graphene sheets. This difference is ex¬ 
plained in terms of the orbital delocalization and barrier vari¬ 
ation for the two geometries. 

The paper is organized as follows. An overview of the the¬ 
ory and the computational approach are given in Section 11. In 
Section III, we discuss the calculation of transport properties 
of graphene/graphene side-contact junctions and present the 
results. Conclusions are in Section IV. A detailed derivation 
of the method is provided in Appendices. 


II. MODEL AND METHODS 
A. Overview of the method 

We consider a system as sketched in FIG. [T] consisting of 
a central scattering region connected to left and right semi¬ 
infinite electrodes. Mapping this to the graphene side-contact 
junction, the two electrodes are two semi-infinite graphene 
sheets with bulk potentials, and the scattering region con¬ 
tains the overlapping region plus a few extra layers out¬ 
side the overlapping region on both sides to ensure conver¬ 
gence. The plane-wave multiple scattering theory method 
produces either complex band structure for each of the bulk 
electrodes or the transmission coefficients for the electrode- 
device-electrode assembly. The calculation is divided into two 
stages. In the first stage, the complex band structures of both 
electrodes are calculated. Even though the basic formalism 
for this stage is equivalent between the new method and the 
Choi-Ihm method, the use of the multiple scattering theory 
speeds up the new method by about 17% on a single proces¬ 
sor for the same supercell size and allows it to be more effi¬ 
ciently parallelized. A further speedup is achieved by recog¬ 


nizing that the complex band structures of the electrodes can 
be computed using primitive cells only, and be folded into a 
larger supercell that matches the transverse dimension (per¬ 
pendicular to the current) of the scattering region at the end 
of this stage (Appendix B). In the second stage, a set of linear 
equations are solved to obtain the transmission coefficients. 
While the final step of solving the linear system is identical to 
the Choi-Ihm method, most of the computational time is spent 
in the steps needed to set up the final equations. Again imple¬ 
menting these steps using scattering matrices makes the new 
method significantly more efficient for parallel computation. 

To set up the equations, we first apply in each region the 
two-dimensional periodic boundary conditions in the trans¬ 
verse (x, y) directions, which are perpendicular to the di¬ 
rection of transport. Systems that are not periodic along 
the transverse directions can be approximated by sufficiently 
large two-dimensional supercells with periodic boundary con¬ 
ditions. The Kohn-Sham equation 36 within the framework of 
the plane-wave pseudopotential method is, 
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where Vi oc ( r) is the total screened local potential that in¬ 
cludes the local part of ionic pseudopotential, electrostatic 
and exchange-correlation potential due to valence electrons; 
W“ ( r) are a set of projector functions associated with atom 
a at position r“, and form the nonlocal part of the pseudopo¬ 
tential with the coefficients We use subscript _L to 

indicate vectors in the xy plane, k and R are wave vectors 
and lattice vectors in the xy plane, respectively. 

Following Choi and lhrrP3 we divide each scattering or 
electrode region into slices perpendicular to the transport di¬ 
rection (z). The slices are sufficiently thin so that within each 
slice the local potential can be treated as independent of z. If 
one such region starting from zo — 0 and ending at zn - d is 
divided into N slices, then Eq. 0 can be rewritten for each 
slice (labeled by superscript p) as, 

EVi r) = - £ v V(r) + V£ c (rJ^(r) 
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The solutions of the above inhomogeneous differential equa¬ 
tion for each slice are to be connected together using the 
boundary conditions that the wave functions and their deriva¬ 
tives are continuous to yield the total wave function 0-(r) of 
the entire region in Eq. Q. The first step is to find a set of 
basis functions, 

ltf(r) = <( r ± )e ±ikP ^ (4) 

satisfying a homogeneous equation obtained by setting all 
Calm’s in Eq. 0 to zero, where n labels different homo¬ 
geneous solutions and goes up to Nzd, the cutoff number 
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FIG. 1. (Color online) Schematic (a) top and (b) side views of a graphene side-contact junction. The two graphene layers are connected by 
one sheet overlapping another to form a bilayer boundary region. The system is periodic along the x direction (six supercells are drawn in the 
top view) and has a large vacuum gap (more than 15A) along the y direction. Along the z direction, the system contains the overlapping region 
and several electrode (graphene) buffer “layers” to ensure the convergence of the potential at the boundary to that in bulk electrodes. The unit 
cell contains one edge carbon atom on each graphene sheet, denoted as Cl (top layer) and C2 (bottom layer), respectively. 


of the reciprocal lattice vectors G ± in the xy plane, k p = 
y/2m(E - E p )/Ti is the n' h wave vector along z direction in the 
p' h slice, E is the total incident energy and E p n is the energy 
eigenvalue of the following equation, 
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where cp p (G±) is the Fourier transform of e/> / , ,’(r_). In addition 
to the homogeneous solution basis set, we also need a partic¬ 
ular solution basis set as described below. Starting from Eq. 
(|3]i we set only one of s to one and all others to zero. 
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which yields a particular solution, i// P lm ( r), for each (a, l, m). 
Both the homogeneous and the inhomogeneous solution basis 
sets satisfy the Bloch boundary condition within the xy plane, 
^n(.aiyJj + RJ = e lk±R± i// n (aim)( r). The general solution of 
Eq. ([3b is written as a linear combination of if/ p and 

<(r) = + £ B p nk f n ( rj e -^-^ 

n n 
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where &,■ labels the wavevector of the Bloch and evanescent 
states when computing complex band structure, and labels the 
wavevector of the incident waves of the entire scattering re¬ 
gion when computing the transmission coefficients. 


The coefficients A p , and B p k , which depend on k t , are the 
only unknowns in the wave function and are determined by 
matching the wave functions between adjacent slices, usually 
in a transfer matrix 39 formulation. However, a transfer ma¬ 
trix approach is usually numerically unstable because of the 
appearance of the exponential factors in Eq. (|7|. For k p con¬ 
taining nonzero imaginary parts (corresponding to evanescent 
states), these terms are exponentially decaying and growing 
waves. Therefore, some transfer matrix elements, as given by 
Eq. ( A17| >, will grow exponentially while some others will de¬ 
cay exponentially, creating a numerically unstable system. As 
iteration proceeds, information for the decaying modes will be 
lost, causing the numerical solution to diverge. To avoid this 
numerical instability, we separate the waves into forward and 
backward waves. The forward waves are those that propagate 
or decay in the positive z direction, the backward waves are 
those that propagate or decay in the negative z direction. By 
iterating both types of waves along their exponentially decay¬ 
ing directions we can avoid the numerical instabili ty, a prac- 
tice already adopted in a number of previous studies®!!! This 
is accomplished by rearranging the boundary conditions be¬ 
tween the slices in terms of incident waves of slice p + 1 with 
coefficients A p and B p+1 , the outgoing waves of slice p + 1 
with coefficients A p+{ and B p using the scattering matrix^® 
S (p, p +1) [see Eq. ( A18|>] that couples them, 


A p+1 
B p 


Sn(p,p+ 1) Sn{p,p+ 1) 

Sn(p,p + 1 ) S 2 2(p,p + 1 ) 

h a (p,p+l)C 
h b (p,p+l)C 


A p 

B p+1 


( 8 ) 


(A181. 


The coefficient h(p,p + 1) is also defined in Eq. 

Clearly, S(p,p) = /, h a (p,p) = 0 and h b (p,p) = 0. 

Calculating the scattering matrix for each slice is the first 
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step to obtain the total scattering matrix S (l, N) of the en¬ 
tire region. The second step is to stitch the slices together 
by applying a doubling technique similar to that employed 
in the layer-KKR method. In general, one can obtain the 
scattering matrix for a collection of slices m through k by 
combining the scattering matrices for m through n and for n 
through k (m < n < k) using the following multiple scattering 
equation^. 


S 11 (oi,k) = S 11 (o,k)[I - S n(m,n)S2i(n,k)] *5 n(m,n), 

S 12 (m,k) = S 12 (o,k) + Sn(n,k)[I - Sn(m,ri)S 2i(n,&)] -1 X 
S i 2 (m,ri)S 22 (n,k), 

S 21(01, k) = S 21(01,0,) + S 22(01, n)S 2 t(n, k)x 

[I -S 12 (hi, ri)S 21(0, k)]~ l S n(oi,n), 

S 22(01, k) = S 22(01, n)S 22(0, k) + S22(01, n)S2\(n,k)X 

[I - S 12(01, n)S 2 i(n,k)]~ 1 S n(oi,n)S 22(0, k), 
h a (m, k) = h a (n,k) + S n (n, k)[I - S n (m, n)S 2 \(n, /t)] 1 x 
[5 i 2 (m, n)h b (n, k) + h a (m, n)], 
h b (m, k ) = h b (m, n) + S 22(01, n)h b (n, k) + S 22(01, n)x 
S 2 \(n,k)[I - S 12 (01, n)S 2 i(n,k)]~ l X 
[5 12 (in, n)h b (n, k) + h“(m, «)]■ 

( 9 ) 

Each time these equations are applied, the number of slices 
represented by the scattering matrix is doubled, thus the name 
“layer-doubling”. The final scattering matrix S(l,N) is ob¬ 
tained by applying layer-doubling repeatedly until all slices 
are contained within the scattering matrix. 

While the scattering matrix formalism is numerically more 
stable than the transfer matrix formalism, the latter is compu¬ 
tationally faster than the former. To achieve the optimal bal¬ 
ance between speed and numerical stability, we iterate two of 
the transfer matrices In and 1 12 along with the scattering ma¬ 
trices S 21 and S 22 (see Eq. ( | A20[ > ) for the first few doubling 
steps, when the condition number of the transfer matrix In is 
reasonable. As soon as the condition number of In turns bad, 
we switch over to iterate entirely with the scattering matrices 
S 11 , S 12 , 5 21 and S 22 - 

The final inhomogeneous equation is identical in form to 
Eq. ((S), 


a n 


Sn(hN) S 12 ( 1 ,AO ‘ 
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h a (\,N)C 
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Eq. ( fT0| ) provides a direct connection between the total inci¬ 
dent waves of the entire region with coefficients A 1 and B N 
and the total outgoing waves of the entire region with coeffi¬ 
cients A n and B l . 

There are 4 N 2 D + N or b unknown coefficients counting A 
(2N 2 d), B (2N2d) and C (N or b), with N or b being the total num¬ 
ber of nonlocal spheres characterized by W " m . As mentioned 
above, the total wave function of the entire region in Eq. O 
can be obtained by connecting the wave function in each slice 
[see Eq. (j7])], which means that the total wave function can 
be expressed as a function of the unknown coefficients A, B 


and C. Therefore, the definition of C in Eq. G provides N or b 
equations about the unknown coefficients. In addition, Eq. 

(101 gives us another 2fto equations. So a total of 2N2D+N or b 
linear equations are provided by Eqs. (|2ji and ( [T()| >. Addi¬ 
tional 2 N 2 d equations are needed, which are provided by the 
boundary conditions below. Two different types of boundary 
conditions are needed depending on the physical problem un¬ 
der study. For complex band structures, the generalized Bloch 
boundary condition^ are employed along the transport di¬ 
rection (see Appendix B). For the transmission coefficients of 
the scattering region, continuity conditions for the wave func¬ 
tion and its derivative are applied at the interfaces between 
the scattering region and each side of the electrodes (see Ap¬ 
pendix C). 

To compute the complex band structure of an electrode re¬ 
gion, the generalized Bloch conditions along the ’-direction 
[see Eqs. CBT), (|B2|) in Appendix B] need to be applied as the 
boundary conditions. With Eq. (10 1 , we can express A N and 


B l as a function of A 1 , B N and C. In addition, we can denote 
the coefficients C with nonlocal spheres completely fitting the 
electrode region as C( a i m y, and those with nonlocal spheres 
crossing the boundaries of the electrode region as C( a i m ). From 
the definition of C( a i m y in Eq. G- we can express C( a [ m y as 
a function of A N , B 1 , A 1 , B N and C( a i m ). Therefore, Eq. (lOl 
and the definition of C( a i m y in Eq. B give us an expression 
of A n , B l , C( a [ m y as a function of A 1 , B N and C( a i m ). So in the 
definition of C ( a i m ) [Eqs. (|B3j) and ( |B4[ )] and Bloch condition 
Eqs. (Bl 1 and (B2 1 , we can substitute A N , B l , C( a i m y wit h A 1 , 
B n and C (r,im i- As a res ult, the only unknowns in Eqs. (B3( 
and (|B4|i, and (Bl 1 and (B2 1 are A 1 , B N and C u ,i ml . Rearrang¬ 
ing these equations by setting the unknowns A 1 , B N and C(„h nj 
as X, we can obtain the generalized eigenvalue problem that 
takes the form 


P X = e ,fc V k| ‘ a3l 2 • X, 


( 11 ) 


where a 3 is the third lattice vector which is not in the xy plane; 
k = (kj_,fc) is the wave vector. Solving this equation yields 
a set of forward waves {^(r)[, backward waves {i/q)(r)[ and 
also complex band structure k(E). 

For transmission coefficient calculations, we need to match 
the boundary conditions at zo = 0 and zn = d for the wave- 
function and its derivative [see Eqs. ( |C14| i-( |CT6| i in Appendix 
C], Through this process. A 1 and B' c an be expressed as a 
function of A 0 , B° and C [see Eq. (C5 1 ]; A N and B N can be 
expressed as a function of A N+1 , B N1 1 and C [see Eq. (C 61 ], 
where A 0 and //', ,4 :V+1 and /f v+1 are the wave coefficients in 
the 0 ,h and (N + 1 ) ,h slices, as defined in Appendix C. There¬ 
fore, the_unknowns in the definition of C of Eq. 0, and Eqs. 
(C14i-(C16i are 4 :V+1 , B {] and C (noting that A 1 ' and B N+ 1 


specify the incident waves, which are known). Rearranging 
these equations by setting the unknowns A :V+I , B° and C as X, 
we can write a set of linear equations in the matrix form 


M ■ X - D. 


( 12 ) 


The reflection and transmission matrices can be obtained from 
the solution of X in Eq. (12 1 . The details about the dimension 
of matrices P, Q, M, D as well as the elements of X are dis¬ 
cussed in Appendices B and C. 
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Next we compare the efficiency of this method against that 
of the Choi-fhm method, fn the first stage, in which the com¬ 
plex band structures of the electrodes are calculated, there 
are two speedups. First, by using the primitive cells of the 
electrodes, the dimension of Eq. (11 1 for complex band cal¬ 
culations can be reduced by up to several fold for a large 
calculation. This strategy can also be implemented in the 
original Choi-Ihm method, although it is not yet available in 
PWCOND. Second, on a single processor, the new method 
speeds up the complex band calculations by about 17% ac¬ 
cording to benchmark runs. There is an additional speedup 
in a parallel calculation. Since this parallel speedup is the 
same for both stages, we describe it specifically for the sec¬ 
ond stage. The most time consuming parts for either stages 
are the steps to set up Eq. ( pT| ) for the complex band calcu¬ 
lation or Eq. (12 1 for the transmission coefficient calculation, 
not the final step of solving either equation (which is paral¬ 
lelized efficiently but does not consume much computational 
time). Both the new method and the Choi-Ihm method have 
identical Eqs. (|5]» and (|6| used at the beginning of the calcula¬ 
tion, for which nearly perfect parallelization can be achieved. 
The major difference between the two methods is in the steps 
immediately following this calculation. In the multiple scat¬ 
tering method, the scattering matrix for each slice [Eq. © is 
computed, then the doubling technique [Eq. (J9)] is applied to 
obtain the scattering matrix [Eq. < p~0] >] for the entire system. In 
the Choi-Ihm method, the wave functions in individual slices 
are matched across the boundaries in a layer-doubling process 
until the wave function of the entire system is obtained. In 
both methods the bottleneck for parallelization is the layer¬ 
doubling step, which can scale at best as log, L where L is 
the number of processors used. However, while the Choi-Ihm 
method relies entirely on the layer-doubling technique for the 
wave functions, in the multiple scattering method much of the 
computation load is shifted to the extraction of the scattering 
matrices for individual layers, which can be parallelized with 
nearly 100% efficiency. Once the scattering matrices of all 
slices are calculated, evaluating Eq. (|9]) is twice faster than 
matching boundary conditions for the wave functions across 
layer boundaries. Therefore the prefactor of the log 2 L term 
for the new method is half of that for the Choi-Ihm method. 


B. Models and computational details 

We apply the plane-wave-multiple-scattering trans¬ 
port method to study the transport properties of a 
graphene/graphene side contact junction, as illustrated 
in FIG. |T] The edges of the graphene sheets are zigzag 
H-terminated in the supercell, the unit cell contains one 
edge carbon atom for each sheet with a dangling bond 
saturated by the H atom. We denote these edge atoms as 
Cl (top graphene layer) and C2 (bottom graphene layer), 
respectively. Transmission and reflection coefficients of 
the incident electrode Bloch states through the scattering 
region are calculated by applying the scattering boundary 
conditions, Eqs. ( |C14|C16| . The electronic structures of the 
electrodes and the scattering region are calculated separately 


using the QUANTUM ESPRESSO package to obtain the 
self-consistent potentials. Fourteen unit cells of graphene 
outside the overlapping region are added to the scattering 
region on each side to ensure the convergence of the potential 
at the boundaries of the scattering region (see FIG. [TJ. The 
electronic structure calculation for the scattering region is 
obtained by repeating the region in the z direction. To achieve 
translational invariance under an unit lattice vector along the 
z direction, the graphene sheet on the left must coincide with 
the graphene sheet on the right after translation. The best way 
to accomplish this is to tilt the supercell such that the lattice 
vector a 3 of the scattering region is no longer perpendicular 
to the xy plane (see FIG. [TJ)). Such a configuration cannot 
be handled by the existing PWCOND code as we discussed 
earlier. The PWCOND code can only treat this system by 
doubling the size of the supercell to make the repeating lattice 
vector perpendicular to the xy plane, thus incurring an 8-fold 
increase in computation time and 4-fold increase in memory 
requirement. 

We use the Rappe-Rabe-Kaxiras-Joannopoulos ultrasoft 
pseudopotentials® for the C, H atoms with the Perdew- 
Burke-Ernzerhof exchange-correlation functional 47 . Because 
of the zigzag H-terminated edges which have non-zero 
magnetization 4 -! the calculation is spin-polarized. The en¬ 
ergy cutoffs are 60 Ry and 360 Ry for the wave function and 
charge density, respectively. Gaussian smearing with a width 
of 0.05 eV is used for the energy levels. The scattering region 
is periodic along the x direction with a thick vacuum layer of 
more than 15 A in the y direction. A 2.46 A x 20 A rectangu¬ 
lar area is used for the supercell’s xy plane. For self-consistent 
calculations, the z dimension of the supercell is much larger 
because of the added unit cells of electrodes, and also depends 
on the width of the overlapping region. A 20 x 1 x 1 k-space 
mesh is sufficient to sample the Brillouin zone (BZ) for the 
scattering region. 

Graphene has a zero density of states (DOS) at the Dirac 
point. Using it as electrodes requires the transport calculations 
to be performed at an energy away from the Dirac point to en¬ 
sure a finite transmission. In our calculation the incident elec¬ 
tron energy is chosen to be 0.1 eV above the Dirac point. At 
this energy, electrode conduction channels only exist within a 
small volume of the reciprocal space. The reciprocal space is 
discretized as k ± = (m x Ak x , m y Ak v ), where A k x (y) is the mesh 
interval along x(y) direction, and m x(y> runs from N x(} 1 to N^ y> . 
The number of k ± points along the x(y) direction in the recip¬ 
rocal space is N x< - y> = (V 4 ' ’ - N xt -‘ 1 + 1, which defines the num¬ 
ber of conduction channels in the electrode. The transmission 
coefficient T(m x Ak x , m y Ak y ) for each mesh point is computed 
using the method described in Section [II] Because the super¬ 
cell along the y direction is the vacuum layer of y direction in 
the electrode in our model system is sufficiently large that the 
total energy is independent of k y , the transmission coefficients 
are also independent of k y . 

Because of the small number of conduction channels in 
graphene, even if the junction does not contain any scatter¬ 
ing the total transmission is still small, a situation we refer to 
as “electrode-limited” conduction. In order to distinguish the 
effect of junction scattering from that of the electrode-limited 
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conduction, we compute a transmission probability per elec¬ 
trode conduction channel as follows: 

N} N> f 

f =WN~y Z E nm*Ak x ,myAk y ). (13) 

m x =Nf m y =N? 

Eq. ( [T3] l can be reduced to the following form since the trans¬ 
mission coefficients are independent of k y : 

Nf 

T = Jp Z T (m x Ak x ,0). (14) 

m x =N? 

We find that A k x = 4.167 x 10~ 4 (27r/ai) is sufficient to con¬ 
verge f to 1%. Here a\ is the lattice constant of electrode 
in x direction. In particular, if the k mesh in the reciprocal 
space is infinitesimal or continuous, Eqs. ( fi~3j ) and ( [T4| > can be 
rewritten in an integral form: 


1 

1 T (k x , ky)dk x dk y . 

(15) 




1 f 


T = 

— T(k x ,0)dk x . 

(16) 


/l eff AV 



where the effective projected Fermi area of the electrode, de¬ 
noted as £2 ei f, is the area projected to xy plane in k space from 
the effective Fermi volume of the electrode; A* s is the effec¬ 
tive projected Fermi length in x direction of k space for the 
electrode graphene. 

III. RESULTS AND DISCUSSION 

A side contact formed between two graphene sheets has al¬ 
ready been observed experimentally and is found to be more 
resistive than grain boundaries formed between graphene do¬ 
mains within the same sheef®. The conduction through the 
side contact can be greatly improved by increasing the over¬ 
lapping area in the experimental samples. Therefore, the ef¬ 
fect on the transport properties of a side contact junction due 
to the overlapping area as a function of the interlayer distance 
between the two graphene sheets is the focus of this first- 
principles study. We begin with the limiting case that the two 
graphene domains have no overlapping with each other. In 
this case the edges play an important role, where the domi¬ 
nant factor is the orbital hybridization between the states from 
the edges of the two graphene layers, which varies with the 
interlayer distance. This hybridization affects the degree of 
localization of the electron orbitals near the edges, thus have 
a deciding role on the transport properties. 

A. Tunneling without overlapping between two graphene 
layers 

We first examine how the interlayer distance (d ± ) and hori¬ 
zontal distance (d\\) affect the tunneling properties when there 
is no overlap between the two sheets (see FIG.[2ji). The spin- 
polarized calculations show that the edge of each graphene 


layer has a small amount of magnetization and the anti¬ 
ferromagnetic (AFM) bonding state between the two edges 
is the ground state. The carbon atoms on equivalent positions 
from the two edges in each layer, e.g.. Cl and C2, have the 
same magnitude of magnetization with opposite signs. Due to 
this spatially antisymmetric spin configuration, the transmis¬ 
sion coefficient for the spin-down channel is identical to that 
for the spin-up channeP®. 

With fixed c/ ± , the calculated transmision as a function of d\\ 
shows an exponential decay, as plotted in FIG. [2J? for d ± = 0 
and d ± =2.1 A. Such a result is expected for simple tunneling 
through a vacuum barrier with a thickness equal to the sepa¬ 
ration d\\. 

The dependence on the interlayer distance (d ± ) is more 
complex. We set d\\ — 0 which yields the maximum trans¬ 
mission for non-overlapping sheets. In this case, the transmis¬ 
sion coefficient first increases with the distance, reaching its 
maximum at d ± =2.1 A, followed by an exponential decrease 
for large distances (see FIG. [ 2 J:). The distance of 2.1 A co¬ 
incides with the minimum of the total absolute magnetization 
as shown in FIG. & which also plots the total energy as a 
function of d ± . The energy difference between ferromagnetic 
(FM) and AFM states is 5 meV when d ± = 1.5 A. For other 
interlayer distances, we cannot find the FM states in the cal¬ 
culation. The equilibrium is reached when the interlayer dis¬ 
tance is about 3.0 A. The fact that the maximum transmission 
coincides with the minimum of absolute magnetization sug¬ 
gests that there is a “cancellation” of the magnetic moments 
between the two opposing edges due to the majority spin elec¬ 
trons from each edge “leaking” into the minority spin channel 
of the other side. 

To understand the mechanism of this spin leakage, we ex¬ 
amine the degree of the orbital hybridization between the edge 
atoms of the two graphene layers. In FIG. [3] we plot the pro¬ 
jected density of states (PDOS) of edge Cl and C2 atoms at 
several interlayer distances. When the interlayer distance is 
large (e.g. d ± = 4.0 A, FIG. [3^), Cl has a peak below the 
Fermi level for spin-up and a peak above the Fermi level for 
spin-down and C2 is the opposite, which gives them the same 
magnitude of magnetization but with opposite signs. When 
the distance becomes smaller (below around 3.0 A), we ob¬ 
serve a small peak above the Fermi level in the spin-up PDOS 
of Cl (see FIG. [ 3 J 5 ), and correspondingly a small peak above 
the Fermi level in the spin-down PDOS of C2 (not shown 
here). This clearly indicates hybridization between the or¬ 
bitals of Cl and C2. The degree of hybridization between 
the orbitals of the two graphene edges can be estimated from 
the size of the small peak in the spin-up PDOS of Cl. The 
small peak reaches its maximum when the interlayer distance 
is 2.1 A, indicating that the orbital hybridization is strongest 
at this distance, leading to the largest transmission coefficient. 
The size of the small peak is almost the same for d±_ = 1.5 A 
and 3.0 A, consistent with Fig. [2J: which shows that the trans¬ 
mission coefficient at these two interlayer distances is almost 
equal. 

The degree of orbital hybridization as a function of d ± can 
be visualized directly from the spin-up integrated local density 
of states (LDOS), as shown in FIG.[4]'a-d), which is calculated 





7 


(a) 



(b) (c) 





FIG. 2. (Color online) (a) Schematic side view of two layer graphene with horizontal distance d\\ (parallel with graphene plane) and interlayer 
distance d ± (perpendicular to graphene plane). The supercell is tilted so that the left of the bottom layer graphene can match the right of the 
top graphene layer with a periodic boundary condition, (b) Transmission as a function of d\\ with d ± = 0.0 A and 2.1 A. (c) Transmission for 
spin-up channel in the left panel and total energy (£), absolute magnetization (M) in the right panel as a function of d ± with d^ = 0. y axis in 
both b and c is in logarithm scale. 


by integrating the spin-up DOS from 0.05 eV below the inci¬ 
dent energy to 0.05 eV above the incident energy (in compli¬ 
ance with our smearing parameter). These plots show that the 
orbitals around incident energy that carry current are mainly 
the carbon p z orbitals. When r/ ± =4.0 A (FIG.[4|j), there is neg¬ 
ligible LDOS on the top graphene layer (contains Cl) indicat¬ 
ing that there is almost no hybridization between the carbon 
orbitals from different graphene sheets. This explains the ex¬ 
ponentially decaying part of the transmission curve at large in¬ 
terlayer distance. With decreasing interlayer distance, one can 
observe some spin-up LDOS on the top graphene layer due to 
the hybridization between the spin-up orbitals of two carbon 
atoms at the edges of two sheets. Because of the hybridization, 
some electrons with the spin opposite to the magnetization di¬ 
rection appear on the edge, as if they are “leaked” from the 
other edge whose magnetization is in the opposite direction, 
reducing the total absolute magnetization of the junction (see 
FIG.[2j;). The p z orbital in Cl, which has the largest spin leak¬ 
age compared to other carbon atoms on the top graphene layer, 
tilts with different angles at different interlayer distances, as 
shown in FIG. &■ At an angle when the orbital points di¬ 
rectly at the C2 atom on the other edge, the hybridization, 
spin-leakage and degree of the LDOS delocalization in the top 


layer graphene are the largest. This happens at the interlayer 
distance of 2.1 A (see FIG.|4})). 


B. Tunneling between two overlapping graphene layers 

In the case of two graphene layers overlapping each other, 
transmission depends sensitively on how the two layers are 
stacked together, and even whether there is a rotation angle 
between the two layer^ 12 ! There are two typical types of stack¬ 
ing. One, in which every atom on the second layer lies over 
an atom of the first, is called the AA stacking. The other, in 
which half of the atoms in the second layer lie directly over 
the center of a hexagon in the lower sheet and the other half 
over an atom, is called the AB stacking. We first examine 
the AA stacking pattern. We fix the interlayer distance of 
graphene layers to that of bilayer graphene, which is about 3.4 
ApSi. By varying the overlapping area, expressed in the unit of 
graphene primitive unit cell area (So), we can calculate the 
transmission coefficient as a function of the overlapping area, 
which is plotted in FIG. [5ji for AA stacking. The transmis¬ 
sion first increases superlinearly with the overlapping area, as 
it varies from two to eight graphene primitive cells. Then the 















8 


(a) 


T^I.O 

c 

'a. 

« 0.5 


> 

<1) 

CO 

® 

re 


0 


<n o.5 

(/) 

o 

Q 1.0 


-1.5 -1 -0.5 0 0.5 1 

Energy (eV) 


1.5 


(b) 


=J 

J2 


< 

w 

o 


Q 



Energy (eV) 



FIG. 3. (Color online) (a) Projected density of states (PDOS) for the p orbitals of the edge carbon atoms. Cl and C2, for d ± = 4.0 A. At this 
distance, there is little evidence of coupling between the two graphene sheets, (b) Spin-up PDOS for the p orbital of Cl at different interlayer 
distances. Here we use a smearing of cr = 0.005 Ry, and the Fermi energy is at 0 eV. 



FIG. 4. (Color online) Isosurfaces (at 0.0002 states/A 3 , in red 
color) of integrated spin-up local density of states from 0.05 eV be¬ 
low the incident energy to 0.05 eV above the incident energy for 
two graphene layers with interlayer distance (a) 1.5 A, (b) 2.1 A, (c) 
3.0 A, (d) 4.0 A. (e) The tilting angle of the p z orbital in Cl of the top 
sheet as a function of the interlayer distance. Yellow balls are carbon 
atoms and blue balls are hydrogen atoms. 


increase slows for larger overlapping areas until the transmis¬ 
sion coefficient appears to converge when the overlapping area 
exceeds ten graphene primitive cells. For large overlapping 
areas, electrode-limited conduction is reached, evident from 
the transmission per conduction channel close to unity, a sit¬ 
uation for which Eq. (13 i is designed to uncover. When the 
overlapping area is small, it presents a constriction to the cur¬ 
rent in the junction. The transmission per conduction channel 


of electrode is expected to be much smaller than unity in this 
case, and it is confirmed by our calculation. Small transmis¬ 
sion due to the constriction at the junction will be referred to 
as junction-limited transport. 

For AB stacked junctions, the transmission as a function 
of overlapping area is shown in FIG. [5j>. The transmission 
for AB stacking is about an order of magnitude smaller than 
that for AA stacking, as shown in FIG. [5]3. Consequently 
the transmission per conduction channel is much smaller 
than unity, placing AB stacked junctions within the junction- 
limited regime. The linear dependence of the transmission on 
the overlapping area depicted in FIG. [ 5 J 5 indicates a simple 
scaling of the tunneling current with area. This may also be a 
consequence of weak coupling between the two layers, as ev¬ 
ident from the much smaller the total integrated LDOS of the 
carbon p- orbitals at the incident energy that carry the current 
for AB stacking than for AA stacking, as shown in FIG. [6] 

To compare and relate the transport properties between AB 
and AA stacked graphene junctions, we examine the effec¬ 
tive potential that electrons experience within the overlapping 
area between two graphene layers, plotted for the region out¬ 
side the core radius, where the pseudopotential is equal to the 
all-electron potential. FIG. [ 7 ] shows the isosurface of the ef¬ 
fective potential at the incident energy (red color). The region 
of interest is between the two layers. The volume is divided 
by the red surface into regions with higher potential than the 
incident energy, which are enclosed by the red surface and ap¬ 
pear as solid red volumes, and regions with lower potential, 
which are the white colored regions between the two sheets 
(those outside both layers also have higher potentials). In the 
case of AA stacking, regions with higher potentials between 
two layers are contained within isolated pockets, and regions 
with lower potentials form connected paths through the entire 
interlayer volume. For AB stacking, the white color low po¬ 
tential regions between the two layers are blocked off by the 
red color high potential regions. This difference in the topog¬ 
raphy of the potential can cause large difference in the trans- 
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FIG. 5. (Color online) Transmission as a function of overlapping area for (a) AA stacking and (b) AB stacking. 
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FIG. 6. (Color online) Isosurfaces (at 0.001 states/A 3 , in red color) 
of integrated total local density of states (summed over both spins) 
from 0.05 eV below the incident energy to 0.05 eV above the incident 
energy for (a) AB stacking with overlapping area 75o, (b) AA stack¬ 
ing with overlapping area 85 o- 5o is the area of graphene primitive 
unit cell. Yellow balls are carbon atoms and blue balls are hydrogen 
atoms. 



(a) 


FIG. 8. (Color online) Transmission as a function of interlayer dis¬ 
tance for both AA stacking with overlapping area 85 o (red) and AB 
stacking with overlapping area 75 o (blue). 5o is the area of graphene 
primitive unit cell, y axis is in logarithm scale. 




FIG. 7. (Color online) Isosurface of the effective potential that elec¬ 
trons experience at the incident energy for (a) AB stacking with over¬ 
lapping area 75o, (b) AA stacking with overlapping area 85o- So is 
the area of graphene primitive unit cell. The isosurface of the effec¬ 
tive potential is in red. Yellow balls are carbon atoms and blue balls 
are hydrogen atoms. 


mission of the electron wave function at the incident energy, 
leading to both differences in the conductance as well as wave 
function hybridization. 

We also plot in FIG.[8]the transmission coefficient as a func¬ 
tion of the interlayer distance between the two layers for both 


AA and AB stacking. When the interlayer distance is larger 
than 4.5 A, both AA and AB stacking reach the vacuum tun¬ 
neling regime giving us the same decaying rate. But the AB 
stacking reach the vacuum tunneling regime at a shorter inter¬ 
layer distance than AA stacking, confirming that the interac¬ 
tion between the two layers is much smaller for AB stacking. 

IY. CONCLUSION 

By implementing the multiple scattering theory within the 
plane-wave basis, we have improved over previous plane- 
wave based transport calculation in terms of both speed 
and parallel efficiency. We apply this method to study a 
graphene/graphene side contact junction system where the 
contact is formed by stacking two graphene layers through 
van der Waals interaction. The transmission through such a 
junction is closely related to spin leakage between the two 
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graphene edges, a consequence of orbital hybridization be¬ 
tween the carbon atoms across the layers which leads to 
the delocalization of the DOS. When the overlapping area is 
large, stacking pattern becomes an important factor in decid¬ 
ing transport properties across the layers. The transmission 
coefficients for AB stacking is one order of magnitude smaller 
than that for AA stacking, primarily due to the larger volume 
of the blocking potential within the overlapping region for AB 
stacking. 
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Appendix A: Layer multiple scattering theory for nonlocal pseudopotentials 

The general solution of Eq. 0 can be expressed as 

*A(r) = ^ ^ ^n^n(X) + ^ ] C almlffaim (O. 


aim 


where and i// a [ m ( r) are the solutions of the following homogeneous and inhomogeneous equations, respectively, 

n 2 


2 m 
2m 


V“(A„(r) + Vi oc (r)i// n (r) = Eif/ n (r), 

V 2 i/Wr) + Vi oc (r)if/„, m (r) + ^ e' k± ' R± W^,(r - r“ - RJ = £i/Wr). 


(Al) 

(A2) 

(A3) 


We divide the system along the z-direction into N slices. If the slices are sufficiently thin then within each slice V/ oc can be 
approximated as independent of z. Applying Fourier transformation for the wavefunctions in the p ,h slice, we have. 


Vni r) = K(rJe ±ikP ^ = £ tf(G 

G ± 

= Y = Y Z <(Go e ,ckj+Gjr ‘/u«’ 


(A4) 

(A5) 


Gx j 


where G ± ) obeys the following eigenvalue equation. 


E P <P P (G ± ) = ^|k ± + G ± U p (GJ + 2 - Gl)^(Gl), 


(A6) 


and fiaiJd is s iven b y- 


fUJz) = GJYe-^^ f '' dz'g p (z - z')W? m 

c .. Jzp-i 


(k ± + Gj_,z' - tT), 




d 2 r ± W°,(r)e- i(k±+G ^ 


(A7) 

(A8) 


Here Qzd is the cross-sectional area of the two-dimensional supercell, k p = y2m(£’ - E p )/h is the z-component of the wave 
vector in the p' h slice, V P oc (G±) is the Fourier transform of the local potential in the p ,h slice, g P (z) = exp (ik p ■ z)/2ik p when z > 0 
and g p (z) = exp (~ik p ■ z)/2ik p when z < 0. Calculations of y„(Gj_), <p p (G) and f p , (z) are similar to the Choi-Ihm method 22 ! 

J J J “ J j,aim 

The total wavefunction in the p' h slice in Eq. (|7]» is obtained from, 




(A9) 


Galm,ki — 


^ ^ D/„x,vv ^ 1 I dz I 

MV p= 1 dz p -\ J 


dV[W^ v (r-r“)]*<.(r). 


(A10) 
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Matching the boundary conditions for the wavefunction and its derivative between any two adjacent slices, we obtain the 
recurrence relation for the expansion coefficients A P k and B P nk , 

i 

1 


A p - 
nk < 2 J£ 


aim 


B p = — 
nki 2 kl 


2 exp(-ik p+1 Az) (k p + k p+1 ) 2/ M p + l )*M p 
j Is 

+ J] B p l l exp(ik p+] Az) (k p - k p+l ) ^JM p ^rM p s J + 2 C alm , ki H a nalm (p,p + 1), 

, ( 

2 expf-^f'Az) (k p - k p+l ) Yj M p'T M ns 
j Is 

+ 2 Kk exp( ik p+1 Az) (k p + k p+l ) } + 2 C alm , ki H b nalm (p,p + 1), 


(All) 


aim 


where H a {p,p + 1), H b (p,p + 1) and M p s are. 


(A12) 

(A13) 
(A14) 
(A15) 

(A16) 

where A p , B p and C are matrices with elements \A P nk ), [B P k } and {C*;,}, respectively; I(p, p + 1) is the transfer matrix between 
adjacent slices p and p + 1 in the following form. 


k*jp>p +d - ^ 2 (M r« - cja 

" j S 

" j S 

M ns = ^ J dr±l</> P (r±)r ■ exp[;(k ± + G ± J • rj, 

and G_ s denotes the s ,h reciprocal lattice vector in {Gj_} and Az the thickness of a single slice along the z-direction. 
The above recurrence relations can be rearranged to form an inhomogeneous equation. 


A p ' 


hi(p,P + 1) In(p,p + 1) 

' A p+1 ' 


H a (p,p+ 1)C 

B p 


hi(p,P + 1) hi(p,P + 1) 

B p+1 

+ 

H b (p,p + 1)C 


Inj(P,P + D=^2 (M r>* M - 


exp(-ik p+l Az)(k p + k p+1 ) exp(ik p+1 Az)(k p - k p+l ) 
exp(-ik p+1 Az)(k p - k p+1 ) exp(ik p+l Az)(k p + k p+1 ) 


(A17) 


The wave function is separated into forward and backward waves with the help of the scattering matrix, whose elements are 
calculated from the transfer matrix, 

S n(p,P + 1) = [hi(p,P + l)] -1 , 

S 1 2 (P,P + 1) = -Un(p,P + l)T l h 2 (p,P + 1), 

Sn(p,p + 1) = hi(p,p + l)S n (p,p + 1), 

Si2(p,p + 1) = hi(p,p + l) + h\(p,p + 1)5 \ 2 (p,p + 1), 
h a (p,p + 1) = ~[I n (p,P + 1 )r l H a (p,p + 1), 
h b (p,p + 1) = H h (p, p + 1)+I 2 \(P,P + 1 )h a (p,p + 1). 


(A18) 


For 1 < m < n < k < N , we obtain the following recurrence relations, 

S ii (m,k) = 5ii(n,k)[/ - Sn(m,ri)S 2 i(n,&)] _1 Sn(m,n), 

S 12 (m,k) = S n(n,k) + S\\{n,k)[I - S i 2 (m,n)S 2 i(n,£)] _1 S n(m,ri)S 22 (n, k), 

S 2 i(m,&) = S 2 i(m,ri) + S 22 (m,n)S 2 \(n,k)[l - Suim,n)S 2 i(n,k)]~ l S n(m,n), 

S2iim,k) = S22{m,n)S22(n,k) + S 2 2(m,n)S 2 \(n,k)[l - S n{m,n)S2\{n,k)]~ l S \2(m,n)S2i(n,k), (A 19 ) 

h a (m, k) = h a {n , k) + S n (n, k)[I - S n (m, n)S 2 i(n, fc )] -1 [ 5 12 (m, n)h b (n,k) + h a {m ,n)], 
h b (m,k) = h b {m,n) + S22(/n,ri)h b (n,k) + 5 22 0 «, n)S2\(n,k)[I - S i 2 (m,n) 5 2 i(n, fc)] _1 X 
[ 5 1 2 (m, n)h b (n , k) + h a (m, n)]. 
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Each iteration of the recurrence relations results in the doubling of the number of slices represented by the the scattering matrix. 
This is called the “doubling technique”. At the end of the iteration, we obtain the scattering matrix that represent the entire 
region [see Eq. ©]• To further improve speed, we use a set of recurrence relations based on a mixture of transfer matrices In, 
I \2 and scattering matrices .S’ 21 , S 22 which is faster but numerically less stable. 


I\\(m,k) = [I]\{m,n) + I\ 2 (m,n)S 2 \(n, k)]I\ \ (n,k), 

I\ 2 (m,k) = [7j 1 (m, ri) + I l2 (m,n)S 2 i(n, k)]I\ 2 (n,k) + I\ 2 (m,n)S 22 (n,k), 

S 2 i(m,k) = S 2 \{m,ri) + S 22 (m,n)S 21 (n,k)[I n (m,n) + I ]2 (m,n)S 2 \(n,k)]~ ] , 

S 22 (m,k ) = S 22 (m,n)S 22 (n,k) - S 22 (m, n)S 21 (n, k)[I\\ (m, n) + f\ 2 (m, n)S 21 (n, k)]~' I \ 2 (m, n)S 22 (n, k), (A20) 

H a (m,k ) = [/] 1 (m, n) + I\ 2 (m,n)S 2 \(n, k)]H a (n, k) + I\ 2 (m,n)h b (n,k) + H a (m,ri), 
h b (m,k ) = h b (m,n) + S 22 (m, ri)h b (n, k) - S 22 (m,ri)S 2 i(n,k)[In(m,ri) + Ii 2 (m,ri)S 2 i(n,k)]~ 1 X 
[I\ 2 (m, ri)h b (n, k) + H a (m, n)]. 


We iterate these equations until the condition number of In exceeds a predetermined criterion. Then we switch to Eq. (A19l 
after computing S n and 5 12 from In and / i2 from Eq. (A18 1 . 


Appendix B: Complex band structure calculation 


If the potential is periodic along the ^-direction, the Bloch conditions can be applied as the boundary conditions, 

M + a 3x,z = d) = e ik ^ +ikd i/, ki (.r ± ,z = 0), 

difr ki (r ± +a 3± ,z = d) = j kra , l+ikd di// k ,(r ± ,z = 0) 
dz dz 


(Bl) 

(B2) 


where a 3 ± is the component of a 3 in the xy plane and d is the z-component of 33 . In a transport problem, the wave functions 
(of either electrode) do not extend to infinity in all directions - they match to boundary conditions at the interfaces between the 
electrodes and the scattering region. Thus the requirement that the wave vector k is real is no longer necessary. Solutions with 
complex k (evanescent waves) are now allowed, changing the above boundary conditions to the generalized Bloch conditionfSl 
States with real k's are the propagating (Bloch) states and those with complex k 's are the evanescent states. 

A complication of imposing the boundary conditions, Eqs. © and ( |B2| i, on a nonlocal pseudopotential, is that one must 
account for the nonlocal spheres that cross one boundary plane and are thus folded to the boundary plane on the other side of 


the supercell by the Bloch boundary condition. This requires Eq.( A10 1 to be rewritten in the following form when the nonlocal 


spheres (characterized by WV ) cross the left boundary of the unit cell at z = 0 , 


C, 


alm,ki 


= y.D 


lm,uv 


M 


d 2 r. 


TO r-Ol^W + e 


-ikd 


w 


d 1 r. 


[W“,(r-T a -dz)]*ik ki (r) 


(B3) 


and for those crossing the right boundary of the unit cell at z = d, 


C a , m , ki =Yj D >”^ f 0 dz f d 2 rAW: v (r - r a )r^k,{r) + e tkd ^ dz J d 2 r ± [W? v (r - r“ + differ) 


(B4) 


The total number of those spheres is N cros i + N c , 
(right) boundary of the unit cell. 


with N cros i(N crosr ) being the total number of nonlocal spheres crossing left 


The unknowns in Eqs. (IOi, (AlOi, (Bli and (B2i are A 1 , A N , B l , B N , C u ,i m y (nonlocal spheres completely fitting the 
electrode region) and C( a i m ) (nonlocal spheres crossing the boundaries of the electrode region). Some of these unknowns, A N , 
B 1 , and C( a i m y , can be eliminated by expressing them in terms of A 1 , B N and C ( ,,/ m) . The latter are collected as a single vector X , 


X = 


A 1 , 

nk, 

b n . 

nk, 

(1 alm),kj 


c, 


(B5) 


Then the remaining equations from Eqs. (fl0|, (|A10|), (|B1|) and (|B2)» are combined into a generalized eigenvalue problem, Eq. 
CD- The dimension of X is (2 N 2 d + N cr0 si + A, T „ ST )x(number of incident waves). P and Q in Eq. CD are (2 N 2D + N cros i + 
Ncrosr) X (2 N 2D + N cros i + N crosr ) matrices. 

Next, we describe how to fold the complex bands of each electrode from the bigger first Brillouin zone of the two-dimensional 
primitive cell, into the smaller Brillouin zone corresponding to the larger two-dimensional supercell of the junction system. Both 
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the two-dimensional primitive cell and its corresponding supercell is in the transverse plane formed by the first two lattice vectors. 
The third lattice vector a 3 is the same for both the primitive cell and the supercell. Let bj, b 2 be the reciprocal lattice vectors of 
the electrode primitive cell in the plane perpendicular to aa and B |, IL the vectors for the supercell in the same plane. We know 
that the supercell’s first two lattice vectors in real space can be expressed as a linear combination of the first two lattice vectors 
of the primitive cell with coefficients being integers, respectively. The inner product of the lattice vector in real space and its 
corresponding reciprocal lattice vector is a constant, In. Thus we can write b| = cnBj + C 12 B 2 and b 2 = C 21 B 1 + C 22 B 2 , where 
cn, C 12 , C 21 and C 22 are integers. We start from the electrode complex band structure described within the larger Brillouin zone 
of the primitive cell. For each k in the first brillouin zone, the Bloch condition along the third lattice vector 33 is, 

tAkfr + a 3 ) = e' k ' a -’i/'k(r), (B 6 ) 

where k = (k ± , /c-j. The complex band structure calculation can give us a series of eigenvalues k- and eigenstates ikk [see 
Eq. {U}]. The Fourier transform of the electrode potential contains only vectors of the reciprocal lattice of the primitive cell, 
m\>i + nb 2 , where m and n are integers. Thus each k in the first brillouin zone of the reciprocal space of the primitive cell 
couples only to k + ;;/b| + «b 2 ±. Then the eigenstate i/// f is superposition of plane waves containing only the wave vector k 
and wave vectors differing from k ± by mbij_ + ubi, 

0 L = EE @mnk\k± + mbi ± + nb 2± ), (B7) 

m n 

where a m „k are the expansion coefficients; |k ± + wbu. + nb 2 ±) is the plane wave basis. 

To fold k = (k x , k z ) of the primitive cell to K = (K ± , K -) of the supercell, we have the following relationship, 

k = K + m k Bj + «%, (B 8 ) 

where m k and n k are integers. Then the Bloch condition along the third lattice vector a? (supercell and the primitive cell have the 
same 83 ) for the supercell is, 


<AK(r + a 3 ) - iAk+»; 1 Bi+h 1 b 2 ( 




Therefore, the complex band at k in the primitive cell will be folded to the complex band at K_ in the supercell. The eigenvalue 
K z -k z - m k Biz ~ // B2/ and the corresponding eigenstate can be expanded in the new plane wave basis |K ± + MBj ± + Nli 2 ±), 
M and N are integers, as: 

^ = X X a ’ nnk '^ A - + mbl± + ” b2 - L ^ = XI XI + ( mk + mc " + nC21 )BlJ - + C+ WC12 + nc 2 2)B 2 ±). (BIO) 


Appendix C: Reflection and transmission calculation 

For transmission coefficient calculations, we need to match the boundary conditions at zo = 0 and zn = d for the wavefunction 
and its derivative, which can be viewed as adding two more slices. We denote the two added slices as the 0 lh and the (N + 1 ) ,h 
slices with wave coefficients A 0 and B l \ A :V+1 and B N+1 . Then we have 


A 1 - 1 
A ik ~ 2 


X4,X«) f 


U'efl 


■E fi ?».E'"l>’ 


k'.eL 


i | ^(Gx,„z = 0) 

^k'iG ± , s ,z = 0) - — exp(ik:Az) - 1 --- 

' kj J az 

i , dfa(G ±iS ,z = 0) 

^(Gi,„z = 0) - -j- exp(ikjAz)— 1 


dz 


(Cl) 




k'.eR 


^(Gi, s ,t = 0) + — exp(-ikjAz) 


dif/k;(G A _,j,z = 0) 


+ X B XX ( < ) 

k'eL s 

+ X [-flalJ °) exp(-fk' Az)], 


^(Gi,„z = 0) + — exp(-ikjAz) 


dz 

#a;(G j _, j ,z = 0) 
dz 


aim 


(C2) 
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A 


N 
jki 


z<?z«> 


ik'eR 


lJ/k’(G^ s ,Z = d ) ■ 


i d^*j(Gj _, s ,z = d) 
kj dz 


+ JX +1 


U.kt 


Z<"5 


fe(Gj.,„z = <7) ■ 


fc'eL s 

+ Xwf-/t ta w|. 

aim 


i dil/ k ’(G ±tS ,z = d) 
k N dz 

j 


(C3) 


Z^'Z'O 


yt'e« 


^k'(Gj_,s,z = d) + 


i di// k ’ (G ± , s ,z - d) 


dz 




k'.eL 


i/'k' i (G ± , s ,z = d) + 


i di// k '(Gj_ !S ,z = d) 


k N 

j 


dz 


(C4) 


where L and R represent the sets of forward and backward waves. In the two new slices, we expand the wavefunction with 
the generalized Bloch basis (including both propagating waves and evanescent waves). Thus for each incident wave k„ the 
dimensions of A 0 , B° and A ;V+1 , A N+l are N^d + N cros i and N 2 D + N crosr , respectively. Recall that the dimensions of A 1 , B 1 and 
A n , A n are all Nzd- Thus we rearrange the above equations into the following form. 


A 1 


Ini 1-0) 7 12 (1,0) 

A 0 


77 fl (l,0)C 

B 1 


/2i(l-0) 7 22 (1,0) 

B° 

+ 

77 fo (l,0)C 


~ A n ~ 


I U (N,N+1) I n (N,N+ 1) 

A w+1 


H a (N,N + 1)C 

b n 


hi(N,N + 1) l22(N,N+ 1) 

b n+ 1 

+ 

H b (N,N+ 1)C 


where 

1.0) =0, 

</,„( 1.0) - -/ ; U(°)exp(-^ A z), 
Hl alm iN,N + \) = -fl lm {d), 
Hh m (N,N+ 1 ) = 0 , 


(C6) 


(C7) 

(C8) 

(C9) 

(CIO) 


^(10)=^^ 

s 


dif/ k f(G ±s ,z=0 ) 

if/ k '(G^ s ,z = 0) - ^ exp(ik .Az) ' - 

, %(Gi,s,z=0) 

</a;(G ±ii ,z = 0) + jjexpi-ikjAz )— i —gj- 


(Cll) 


7;*;(A7,A7+1) = 

s 


i/'>;(G J _, s ,z = <7) - 
<At;(G x ,.s,z = d) + 


; <9fe (G_L„,,z=d) 
k'J dz 

/ dtk t ’(G ± , s ,z=d) 
k N dz 

J 


In addition, Eq. © can be rewritten as 


/ - 5 12 (1, A0 

A w 


Snd.AO 0 

A 1 


h a (l,N)C 

0 -S22(1,A0_ 

B W 


S 2i (1,1V) -I 

B 1 

+ 

h b (l, N)C 


From Eqns. dC5), (|C6| and (IC13J), we can construct the relationship between A 0 , B° and ,4 iV+l , A 


iJV+l aN+ 1 


/ -S 12 (1,A0‘ 

" / n (Al,lV+ 1) I\2(N,N + 1)' 

A w+1 ' 


0 -S 22 (1,1V) 

I 2 i(N,N + 1) I 22 (N,N+ 1) 

B n+ 1 



Sn(LN) 0 
SidhN) -I 


7n(l,0) 7 12 (1,0) 

A 0 


7 21 (1,0) 722(1,0) _ 

B° 

+ 


H a (l,0) + h a (l,N) - H a (N,N + 1) 
H b ( 1,0) + h b ( 1, N) - H b (N, N + 1) 


x 

C. 


(C12) 


(C13) 


(C14) 
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Because the left electrode and the scattering region share those atoms whose nonlocal spheres lie across the z = 0 plane, we have 
Ncrosi additional equations in the form. 


C, 


(alm),kj 


= 2 >, 

k'eR 


LE 4 O 
( 1 alm),k f ik , i ,ki 


■Zd 

k'.eL 




.0 


LE 

(alm),k , i MJ k'.,k l 


(C15) 


Similarly, for those spheres lying across the z = d plane, we have N crosr additional equations 


C 


(alm),kj — 


Z r RE 4 N+I 

alm),k'jk',ki 

k\eR 


+ 


Z r RE pN +1 

^{aim),k’i k^Jii * 

k' t eL 


(Ci6) 


Here LE (RE) means left (right) electrode. 

The refore fo r each boundary cond iton, altogether we have 2Nid + N or b + N cros i + N crosr equations from Eqns. ( |C14[ i (INid), 
( C15 1 (N cros i), ( CI 6 1 (N crosr ), ( AlO i (N or b)- The number of unknowns for |A W+1 ), {B°) and (C) are [Nid + N cros i} + {Nzd + 
Ncrosr) + \Norb)- Note that {A 0 ) and \B N+1 ) provide the boundary conditions that specify the incident wave, e.g. A° = I, B N+l = 0 
for waves incident from the left electrode for which the corresponding transmission and reflection matrices are {A v+I | and {/J 0 }, 
respectively. A 0 = 0, B N+l = I represents waves incident from the right electrode for which the c orresponding transm issi on and 
reflection matrices are {B°) and {,4' V41 ), respectively. As explained in the main text, rearrange Eqs.(C14i, (C151, (C161 and (AlOi, 
we can obtain a set of linear equations listed in Eq. © where A7 is a (2Nod Ncrosl J" Ncrosr J" Norb)X-(2N2D 4 "Ncrosl Ncrosr 4 Norb) 
matrix, D is (2N2d + N cr0 si + N cr0 sr + N or fc)x(number of incident waves), X is also (2N2d + N cr0 si + N cr0 sr + N 0 ,i)x(number of 
incident waves) and has the following structure, 


X = 


A^ + t 

B° 

C 


(C17) 
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